Quantum Monte Carlo calculations of H 2 dissociation on Si(OOl) 
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We present quantum Monte Carlo calculations for various reaction pathways of H2 with Si(001), 
using large model clusters of the surface. We obtain reaction energies and energy barriers notice- 
ably higher than those from approximate exchange-correlation functionals. In improvement over 
previous studies, our adsorption barriers closely agree with experimental data. For desorption, the 
calculations give barriers for conventional pathways in excess of the presently accepted experimental 
value, and pinpoint the role of coverage effects and desorption from steps. 



The dissociative adsorption of molecular hydrogen on 
the Si(001) surface has become a paradigm in the study 
of adsorption systems. Despite its apparent simplicity, 
more than a decade of extensive experimental and theo- 
retical investigations have not clarified fundamental as- 
pects of the chemical reaction of H2 with this surface. 

Many of the experimental observations are hard to 
reconcile in a unified picture: the sticking probability 
for dissociative adsorption of H2 on the clean surface 
is very small at room temperature suggesting a high 
adsorption barrier; sticking increases dramatically with 
higher surface temperatures On the other hand, 

the nearly thermally distributed kinetic energy of des- 
orbing molecules has lead researchers to the conclusion 
that the molecules have transversed almost no adsorption 
barrier M . Microscopically, these observations were orig- 
inally interpreted in terms of an intra-dimer mechanism, 
where the hydrogen molecule interacts with one single 
dimer of the Si(001) surface |J. However, very recent ex- 
periments have pointed to additional mechanisms involv- 
ing not just a single dimer but nearby dimers The 
existence of highly reactive pathways was first demon- 
strated on steps H or H-precovered surfaces , and ev- 
idence that H2 reacts with two adjacent dimers has also 
now been given for the clean surface ||. In Fig. [IJ the 
intra-dimer (H2*) and two intcr-dimcr pathways at dif- 
ferent coverages (H2 and H4) are schematically shown. 

Theoretically, density-functional theory (DFT) calcu- 
lations performed on intra-dimer and inter-dimer mecha- 
nisms have lead to limited agreement with experiments. 
While correctly predicting the existence of a barrier-less 
H4 inter-dimer reaction path at high coverages (|] , previ- 
ous DFT slab calculations yielded an adsorption barrier 
for the low-coverage H2* and H2 pathways too low to ex- 
plain the small sticking coefficient observed at low tem- 
peratures 0. Desorption barriers from DFT obtained 
within the generalized gradient approximation (GGA) 
were also generally lower than the experimental value 
(2.5 eV, Ref. [0). Additional evidence for a possible in- 



adequacy of DFT-GGA to describe this reaction comes 
from comparison with highly correlated quantum chem- 
istry calculations for small cluster models of the surface: 
For the intra-dimer pathway, these methods obtain val- 
ues for the desorption barrier that are at large variance 
with the DFT-GGA value flW&. 




FIG. 1. Intra-dimer (H2*) and inter-dimer mechanisms at 
low (H2) and high (H4) coverages. The surface configuration 
along the pathway is schematically shown: A circle represents 
a Si atom and a filled circle Si-H. E^ ds , Ef cs and E rxn are the 
adsorption, desorption and reaction energies, respectively. 

In this Letter, we use quantum Monte Carlo (QMC) 
techniques on large cluster models of the surface to accu- 
rately compute the reaction energetics of H2 on Si(001) 
for the intra-dimer and inter-dimer mechanisms, as well 
as adsorption at steps. The reliability of a calculation 
depends both on the level of theory at which electronic 
correlations are treated, as well as on the geometrical 
model used to describe the system. Compared to other 
theoretical approaches, QMC offers the advantage that 
accurate reaction energetics can be computed for rela- 
tively large systems. Our calculations predict reaction 
energies and barriers noticeably higher than those ob- 
tained in DFT with the local density approximation or 
commonly used GGAs. Even though distinctly lower, the 
energetics determined with the B3LYP [[l4| hybrid func- 
tional are the only ones in fair agreement with the QMC 
results. When compared to experiments, the QMC ad- 
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TABLE I. Adsorption (E^ ds ), desorption (Ef cs ) and reaction (E rxn ) energies in eV for H2/Si(001) via the intra-dimer H2* 
mechanism (see Fig. (il), calculated within PW91, B3LYP and QMC. Zero-point energies (ZPE) are not included. 





-gads 


PW91 

Ef s 


Erxn 


gads 


B3LYP 

E dcs 


Erxn 


gads 


QMC 

gdcs 


Erxn 


SigHi2 


0.69 


2.86 


2.17 


0.90 


3.40 


2.50 


1.01±0.06 


3.65±0.06 


2.64±0.06 


Sii5Hi 6 


0.56 


2.65 


2.09 


0.71 


3.20 


2.49 


0.98±0.05 


3.52±0.07 


2.54±0.06 


Si2iH2o 


0.32 


2.31 


1.99 


0.56 


2.90 


2.35 


0.61±0.05 


3.11±0.05 


2.49±0.05 


Si27H24 


0.37 


2.35 


1.98 


0.57 


2.91 


2.33 


0.49±0.13 


3.03±0.13 


2.54±0.13 



sorption barriers for the intra-dimer (H2*) and the low- 
coverage inter-dimer (H2) pathways can explain why the 
sticking coefficient on the clean Si(001) surface is so small 
at low temperatures. Moreover, the QMC desorption 
barriers represent an important input for the interpre- 
tation of experimental results. 

Computational methods. We employ slabs as well as 
clusters to mimic the Si(001) surface. The surface can 
be appropriately modeled with clusters containing only 
a single row of dimers |[15[ : interactions are negligible be- 
tween neighboring dimer rows, while substantial between 
dimers in the same row. Such clusters with one, two, 
three, and four dimers are SigHi2, Sii5Hi6, Si2iH2o, and 
Si27H 2 4. They represent a four layer cut of the Si(001) 
surface with all but the surface atoms terminated with 
hydrogens to passivate dangling bonds. Two separate 
clusters are constructed to model the clean surface and 
the surface after adsorption of a hydrogen molecule; the 
transition state (TS) connecting the the two configura- 
tions is then determined. 

As a starting point, plane- wave pseudopotential calcu- 
lations within DFT were performed to obtain geometries 
and total energies of both the cluster and slab models. 
For details on the construction and the geometry of the 
clusters, see Ref. Q. We optimized all geometries using 
the PW91 functional jl7), which gives a good descrip- 
tion of the structural properties of Si (lattice constant 
error < 1%). Moreover, for cluster models of H2/Si(001), 
the geometries optimized using PW91 and the hybrid 
functional B3LYP were found to be very similar |fl8| , 
and the energetics of the reaction on both sets of ge- 
ometries essentially the same within B3LYP (as shown 
below, B3LYP gives energies very close to our accurate 
QMC results). It is therefore a sound procedure to use 
QMC on geometries obtained from PW91 calculations to 
assess whether a more accurate treatment of electronic 
correlation can change the physical picture. 

We calculate the energetics for various reaction path- 
ways of H 2 on the Si(001) using B3LYP and QMC for 
the cluster geometries, and PW91 for both clusters and 
slab. The many-body wave function used in QMC is of 
the form given in Ref. |L9| (modified to deal with pseudo- 
atoms): 

* = dnD n D n II J ( Ti 3 ' r ' la ' r i a ) ' 

n aij 

Djj and are Slater determinants of single particle or- 



bitals for the up and down electrons, respectively, and 
the orbitals are represented using atomic Gaussian ba- 
sis p(J . The Jastrow factor correlates pairs of electrons i 
and j with each other, and with every nucleus a, and dif- 
ferent Jastrow factors are used to describe the correlation 
with a hydrogen and a silicon atom. The determinantal 
part of the wave function is generated within Hartree- 
Fock or MCSCF, using the quantum chemistry package 
GAMESS pl| . The parameters in the Jastrow factor 
are optimized within QMC using the variance minimiza- 
tion method |2^] and the accuracy of the wave function 
tested at the variational level. The wave function is then 
used in diffusion Monte Carlo (DMC), which produces 
the best energy within the fixed-node approximation (i.e. 
the lowest-energy state with the same nodes as the trial 
wave function) [g3| . All QMC results presented are from 
DMC calculations. 

Intra-dimer mechanism. In the H2* pathway, the hy- 
drogen molecule dissociatively adsorbs on the same Si 
dimer through an asymmetric TS. Multi-reference con- 
figuration interaction (MRCI) calculations Jl3| on SigHi2 
yield adsorption and desorption barriers which are 0.3 
and 0.8 eV higher than the corresponding PW91 values 
listed in Table |. Given the magnitude of the discrepancy, 
it is important to ascertain whether this difference will 
persist in going from the one-dimer model to a more re- 
alistic representation of the surface. QMC gives energies 
in very good agreement with the MRCI results for SigHi2 
and, unlike traditional quantum chemistry methods, can 
be applied to larger models of the surface to study the 
convergence of the reaction energies with cluster size and 
access accurate estimates for the real surface. 

In Table | we list the PW91, B3LYP and QMC ad- 
sorption, desorption and reaction energies via the intra- 
dimer mechanism for the clusters with one, two, three, 
and four surface dimers. The QMC energies are con- 
sistently higher than the PW91 values for all clusters: 
The QMC adsorption, desorption and reaction energies 
are above the PW91 values by about 0.3, 0.8 and 0.5 
eV, respectively. Therefore, the corrections to the PW91 
values due to the incorrect treatment of electronic cor- 
relation are indeed significant, and, interestingly, they 
do not show a noticeable dependence on the size of the 
cluster. The QMC results for Si27H24 have a large error- 
bar but also demonstrate the smooth QMC convergence 
with system size. B3LYP represents a significant im- 
provement upon PW91, giving energies which are much 
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closer to our QMC results, and only lower by about 0.05, 
0.2, and 0.15 eV for Si2iH 2 o- This is in accordance with 
earlier studies using the B3LYP functional for the Si-H 



system 



24 



TABLE II. Adsorption, desorption and reaction energies 
in eV per molecule for H2/Si(001) via the H2 and H4 mech- 
anisms (see Fig. [l]). The Si27H24 model cluster of the surface 
is used. 

Eads -r-idcs T-i 

a -EJa Jirxn 
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Inter-dimer H2 and H4 mechanisms. In the H2 path- 
way, the hydrogen molecule dissociates over two clean 
neighboring Si dimers, yielding a cis configuration with 
two hydrogens bound to two Si atoms at the same side 
of the dimers. For the H4 mechanism, the adsorption oc- 
curs on two neighboring Si dimers which are both already 
covered on the same side with hydrogens (see Fig. [l]). 
The configuration after adsorption consists of two fully 
covered Si dimers. The PW91, B3LYP and QMC ad- 
sorption, desorption and reaction energies for the H2 and 
H4 pathways are presented in Table Q. All calculations 
were performed on the four-dimer cluster Si27H24 and 
the reaction occurs on the two central dimers. For both 
mechanisms, the results show the same trends as for the 
intra-dimer pathway: PW91 significantly underestimates 
reaction energies and barriers, and B3LYP is much closer 
to the QMC results than PW91. For the H4 mechanism, 
PW91 and B3LYP predict no adsorption barrier. Within 
the statistical error, we find this to remain true also in the 
QMC calculation. To explore a possible QMC adsorption 
barrier, we computed the QMC energies for nine geome- 
tries along the path connecting the configuration of four 
hydrogens on two neighboring dimers with the configura- 
tion of two hydrogens in the cis configuration plus a des- 
orbed hydrogen molecule. However, the cis configuration 
does not constitute the energetically lowest adsorption 
geometry of two hydrogens atoms, and thus occurs only 
scarcely. The pairing energy, that is the energy difference 
between the cis configuration and two hydrogens at the 
same Si dimer, is found to be 0.34, 0.42 and 0.54±0.07 
eV in PW91, B3LYP and QMC, respectively [||. 

Discussion and conclusions. In addition to accessing 
the reliability of approximate DFT in describing H-Si 
bonded systems, we can use our accurate QMC ener- 
getics to understand the physics of the interaction of hy- 
drogen with the Si(001) surface. In order to compare 
with the available experimental results for the H2/Si(001) 
system, we need to extrapolate the QMC energies to 



TABLE III. Extrapolated QMC adsorption, desorption 
and reaction energies in eV for the H2*, H2 and H4 mech- 
anisms (see text). a Ref. jL8| (identical AZPE assumed for all 
mechanisms); 6 Ref. §; c Ref. (h}; d Ref. @. 
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the infinite-system limit and include the zero-point en- 
ergy corrections (AZPE). To compute the extrapolated 
QMC results, we start from the PW91 infinite-system 
limit, that is the PW91 slab energy, and add a cor- 
rection for the inaccurate treatment of electronic ex- 
change and correlation estimated from the cluster cal- 
culation aS AE rnrr = E„ 



? QMC 



Jcorr — ^cluster E c Zt 9 cV After including 
the ZPE contribution, the final QMC energy is E QMC 

;iQMC 
J clustcr 



(E 



E™ r ) + AZPE. In Table III, the 



t^PW91 

extrapolated QMC energies are compared with the ex- 
perimental results. 

The PW91 slab reaction energies and barriers are sig- 
nificantly lower than the experimental results and, for 
E^ es and E rxn , the discrepancy is even amplified when 
the negative ZPE's are added. As discussed above, the 
QMC energetics for the clusters are higher than the cor- 
responding PW91 values for all the mechanisms, so the 
correlation corrections are positive and the extrapolated 
QMC results significantly higher than the original PW91 
slab energies. The final reaction energies are compat- 
ible with experiments and the adsorption barriers are 
in much better agreement with the experimental results 
than the PW91 barriers. The larger adsorption barriers 
for the intra-dimer and the low-coverage H2 inter-dimer 
pathways can finally explain why the sticking coefficient 
on the clean Si(001) surface is so small at low temper- 
atures. Moreover, the high reactivity for adsorption via 
the H4 mechanism is corroborated by the QMC results 
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which are consistent with a barrier-less pathway. Pos- 
sibly, the dramatic increase in sticking probability with 
surface temperature is partly due to this pathway: if hy- 
drogen is already present on the surface, an increase in 
surface temperature leads to an increased number of cis 
configurations, thus creating barrier-less adsorption sites 
for sticking. 

Concerning desorption, the picture is more compli- 
cated. Judging from the extrapolated QMC barriers, 
none of the studied mechanisms appears to be com- 
patible with the desorption energy of 2.5±0.1 eV ob- 
served in temperature-programmed desorption experi- 
ments [Io|. However, such experiments could have been 
affected by a small concentration of surface imperfec- 
tions. To exemplify this possibility, we looked at the 
adsorption/desorption of H2 at the Db step edge which 
is modeled as a Si2sH28 cluster constructed from previ- 
ous slab calculations . Experimentally, a small adsorp- 
tion barrier of 0.09 ±0.01 eV is observed while PW91 pre- 
dicts no adsorption barrier, so PW91 desorption and re- 
action energies are the same. Here, we compute only the 
B3LYP and QMC reaction energies for Si2sH28, which 
are equal to 2.54 and 2.81±0.07 eV, respectively. The 
extrapolated QMC value is 2.61±0.07 eV, significantly 
lower than the other desorption barriers, and hence the 
presence of a small number of steps could dominate the 
desorption yield. Our calculations suggest that the hith- 
erto accepted experimental value of the desorption bar- 
rier should be referred to desorption from steps or de- 
fects, rather than to one of the mechanisms discussed 
for ideal Si(001) surfaces. For experiments where contri- 
butions from surface imperfections are carefully avoided 
(e.g. by heating the surface only very locally), the QMC 
calculations predict a slight preference for the H4 mech- 
anism. While the lack of an adsorption barrier along the 
H4 pathway can explain the low kinetic energy of the des- 
orbing hydrogen molecules, one should also keep in mind 
the observed vibrational excitation in desorption p7| , in- 
dicating that some molecules desorb on a pathway with 
an adsorption barrier. We conclude that several mech- 
anisms contribute to desorption, whose relative impor- 
tance depends sensitively on temperature, coverage, and 
surface perfection. 

In this Letter, we presented accurate QMC calcula- 
tions for various pathways of adsorption/desorption of H2 
from Si(001), using large cluster models of the surface. 
For intra-dimer and inter-dimer pathways, we find that 
PW91 significantly underestimates reaction energies and 
barriers, while B3LYP yields energetics in much better 
agreement with the QMC values. Caution should there- 
fore be used when employing PW91 or other GGAs in 
the study of Si-H systems. Finally, the QMC adsorp- 
tion barriers are in close agreement with experimental 
values while the results for desorption call for further 
experimental studies of the activation energy and its de- 
pendence on coverage and surface perfection. 
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